This invention relates generally to the field of seismic prospecting and, more particularly, to imaging of seismic data. Specifically, the invention is a method for controlled-amplitude prestack time migration of seismic data.
In the oil and gas industry, seismic prospecting techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. A seismic prospecting operation consists of three separate stages: data acquisition, data processing, and data interpretation. The success of a seismic prospecting operation is dependent on satisfactory completion of all three stages.
In the data acquisition stage, a seismic source is used to generate a physical impulse known as a xe2x80x9cseismic signalxe2x80x9d that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflected signals (known as xe2x80x9cseismic reflectionsxe2x80x9d) are detected and recorded by an array of seismic receivers located at or near the surface of the earth, in an overlying body of water, or at known depths in boreholes. The seismic energy recorded by each seismic receiver is known as a xe2x80x9cseismic data trace.xe2x80x9d
During the data processing stage, the raw seismic data traces recorded in the data acquisition stage are refined and enhanced using a variety of procedures that depend on the nature of the geologic structure being investigated and on the characteristics of the raw data traces themselves. In general, the purpose of the data processing stage is to produce an image of the subsurface geologic structure from the recorded seismic data for use during the data interpretation stage. The image is developed using theoretical and empirical models of the manner in which the seismic signals are transmitted into the earth, attenuated by the subsurface strata, and reflected from the geologic structures. The quality of the final product of the data processing stage is heavily dependent on the accuracy of the procedures used to process the data.
The purpose of the data interpretation stage is to determine information about the subsurface geology of the earth from the processed seismic data. For example, data interpretation may be used to determine the general geologic structure of a subsurface region, or to locate potential hydrocarbon reservoirs, or to guide the development of an already discovered reservoir. Obviously, the data interpretation stage cannot be successful unless the processed seismic data provide an accurate representation of the subsurface geology.
Typically, some form of seismic migration (also known as xe2x80x9cimagingxe2x80x9d) must be performed during the data processing stage in order to accurately position the subsurface seismic reflectors. The need for seismic migration arises because variable seismic velocities and dipping reflectors cause seismic reflections in unmigrated seismic images to appear at incorrect locations. Seismic migration is an inversion operation in which the seismic reflections are moved or xe2x80x9cmigratedxe2x80x9d to their true subsurface positions.
There are many different seismic migration techniques. Some of these migration techniques are applied after common-midpoint (CMP) stacking of the data traces. (As is well known, CMP stacking is a data processing procedure in which a plurality of seismic data traces having the same source-receiver midpoint but different offsets are summed to form a stacked data trace that simulates a zero-offset data trace for the midpoint in question.) Such xe2x80x9cpoststackxe2x80x9d migration can be done, for example, by integration along diffraction curves (known as xe2x80x9cKirchhoffxe2x80x9d migration), by numerical finite difference or phase-shift downward-continuation of the wavefield, or by equivalent operations in frequency-wavenumber or other domains.
Conversely, other seismic migration techniques are applied before stacking of the seismic data traces. In other words, these xe2x80x9cprestackxe2x80x9d migration techniques are applied to the individual nonzero-offset data traces and the migrated results are then stacked to form the final image. Prestack migration typically produces better images than poststack migration. However, prestack migration is generally much more expensive than poststack migration. Accordingly, the use of prestack migration has typically been limited to situations where poststack migration does not provide an acceptable result, e.g., where the reflectors are steeply dipping.
In some cases, reflector dip can exceed 90 degrees. As is well known in the art, it may be possible to image these xe2x80x9coverturnedxe2x80x9d reflectors using data from seismic xe2x80x9cturning rays.xe2x80x9d Prestack migration techniques must be used in order to obtain an accurate image of overturned reflectors from seismic turning ray data.
There are two general types of prestack migration, prestack time migration and prestack depth migration. The present invention relates to prestack time migration which is used in situations where the subsurface seismic velocity varies in the vertical direction, but can be regarded as approximately constant laterally. Prestack time migration is widely used in the petroleum industry; it is generally considered applicable to most (but not all) prospects.
In practice, prestack time migration is usually approximated by a composite procedure involving successive steps of normal moveout correction (denoted as xe2x80x9cNMOxe2x80x9d), followed by dip moveout correction (denoted as xe2x80x9cDMOxe2x80x9d), followed by zero-offset migration (referred to as xe2x80x9cZOMxe2x80x9d). This sequence of processes is performed on collections of seismic data traces having the same or nearly the same source-receiver separation distance and the same or nearly the same azimuthal orientation (referred to as xe2x80x9ccommon-offset, common-azimuth gathersxe2x80x9d). It is known from experience that the combination of these imaging steps results in a good approximation to prestack time migration in most situations. The approximation tends to be less accurate in situations where there are dipping reflectors in a medium whose velocity varies significantly with depth. In such situations, the inaccuracies manifest themselves in mispositioned images of seismic reflectors, as well as in distorted amplitudes. A principal advantage of the decomposition of prestack time migration into the NMO/DMO/ZOM sequence is that the required algorithms execute rapidly on commonly available digital computers, which results in economical processing.
Kirchhoff prestack migration is the only commonly known method of performing prestack time migration which overcomes the inaccuracies of the NMO/DMO/ZOM decomposition mentioned in the previous paragraph, and which is applicable to common offset gathers. However, Kirchhoff prestack migration, as compared to the NMO/DMO/ZOM decomposition, is generally very expensive and time consuming. For this reason, the seismic data processing industry is continuing its efforts to develop techniques that can be used to produce accurate prestack migrations in an economical manner.
In Kirchhoff prestack migration, the data in a common-offset gather are summed over all source-receiver midpoints (denoted by {right arrow over ("xgr")}). During Kirchhoff summation, the data are time shifted by an amount that depends on the assumed seismic velocity structure in the earth. The time shift is denoted by the symbol tD({right arrow over ("xgr")}, {right arrow over (x)}), which depends on the midpoint location {right arrow over ("xgr")} and the imaged point {right arrow over (x)}=(x,y,z). During summation, the data may be multiplied with a weight w that generally depends on {right arrow over ("xgr")} and {right arrow over (x)}. As explained in Schleicher, J. et al., xe2x80x9c3-D true-amplitude finite-offset migration,xe2x80x9d Geophysics, volume 58, pp. 1112-1126 (1993), the weight w may be chosen in such a way that the migration processing preserves the seismic amplitudes. This is usually considered very desirable, since the amplitudes are widely used in interpreting the processed data.
The Kirchhoff prestack migration process may be expressed mathematically in the form of equation (1):                                           D            M                    ⁡                      (                          x              →                        )                          =                              [                                          ∂                                  ∂                  t                                            ⁢                              ∫                                                                            ⅆ                      2                                        ⁢                                          ξ                      →                                                        ⁢                  w                  ⁢                                      xe2x80x83                                    ⁢                                      (                                                                  η                        →                                            ,                      z                                        )                                    ⁢                                      D                    ⁡                                          (                                                                        t                          +                                                      t                            D                                                                          ,                                                  ξ                          →                                                                    )                                                                                            ]                                t            =            0                                              (        1        )            
where DM({right arrow over (x)}) is the migrated image at the point {right arrow over (x)}, D(t,{right arrow over ("xgr")}) denotes the input data at time t and midpoint {right arrow over ("xgr")}, and w({right arrow over (xcex7)},z) is the weight function discussed above, expressed as a function of {right arrow over (xcex7)} and z (depth). The variable {right arrow over (xcex7)} in equation (1) is the displacement between the input trace midpoint and the surface coordinates of the image point, {right arrow over (xcex7)}={right arrow over (x)}xe2x88x92{right arrow over ("xgr")}. Since, as explained above, offset is held fixed during prestack migration, there is no reference to it in equation (1) or any of the other equations set forth herein.
A sketch of the geometry is included in FIG. 1. A seismic source 10 and a seismic receiver 12 are located on the surface of the earth 14. The midpoint between source 10 and receiver 12 is indicated by {right arrow over ("xgr")}. Incident seismic ray 16 encounters subsurface reflector 18 at reflection point {right arrow over (x)}, and the resulting reflected seismic ray 20 is received by seismic receiver 12. The seismic signal traveltime from source 10 to reflection point {right arrow over (x)} is t1. The seismic signal traveltime from reflection point {right arrow over (x)} to receiver 12 is t2. The total seismic signal traveltime is denoted by tD=t1+t2. As noted above, the horizontal displacement between midpoint {right arrow over ("xgr")} and the surface (i.e., x and y) coordinates of reflection point {right arrow over (x)} is represented by {right arrow over (xcex7)}.
In equation (1), the initial partial derivative indicates that the data should be differentiated with respect to time in order to preserve the shape of the seismic pulse during processing. All quantities in equation (1) are taken to be at fixed source-receiver offset. As more fully discussed below, a difficulty with implementation of equation (1) is that the summation over {right arrow over ("xgr")} generally requires time-consuming computation.
Equation (1) operates in the space-time ({right arrow over ("xgr")}xe2x88x92t) domain. For the special situation of zero-offset migration, it is well known that equation (1) can be implemented more efficiently by operating in the frequency-wavenumber (xcfx89xe2x88x92{right arrow over (k)}) domain. This is described, for example, in Kim et al., xe2x80x9cRecursive wavenumber-frequency migration,xe2x80x9d Geophysics, volume 54, pp. 319-329 (1989).
In Dubrulle, xe2x80x9cNumerical methods for the migration of constant-offset sections in homogeneous and horizontally layered media,xe2x80x9d Geophysics, volume 48, pp. 1195-1203 (1983), the author hypothesized a possible approach to xcfx89xe2x88x92{right arrow over (k)} migration at finite offset, and demonstrated his method applied to a simple example for 2-D data. However, Dubrulle did not discuss how to implement his method for 3-D data. In addition, Dubrulle did not address the issue of how to implement xcfx89xe2x88x92{right arrow over (k)} migration to preserve seismic amplitudes, nor did he address the issue of how to handle seismic turning rays. The present invention overcomes these limitations.
A recent example of a prestack time migration method which operates in the frequency-wavenumber domain and has similar accuracy to Kirchhoff prestack migration is described in Etgen, xe2x80x9cV(Z) F-K migration of common-offset common-azimuth data volumes, part I: theory,xe2x80x9d 68th Annual International Meeting of the Society of Exploration Geophysicists, Expanded Abstracts (Sep. 13-18, 1998). According to Etgen, accurate prestack time migration can be performed using multiplication in the xcfx89xe2x88x92{right arrow over (k)} domain of a migration operator with the data volume.
The method proposed by Etgen is shown pictorially in FIG. 2. The unmigrated common-offset data traces D(t,{right arrow over ("xgr")}) are 3-D Fourier transformed from the space-time ({right arrow over ("xgr")}xe2x88x92t) domain to the frequency-wavenumber (xcfx89xe2x88x92{right arrow over (k)}) domain resulting in transformed data {tilde over (D)}(xcfx89,{right arrow over (k)}). In FIG. 2, this is represented by the transformation from data cube 22 to data cube 24. The migration operator M(xcfx89,{right arrow over (xcex7)},z) is first computed in the xcfx89xe2x88x92{right arrow over (xcex7)} domain for every depth z1, z2, . . . zN. In FIG. 2, the migration operator is represented by cubes 26, 28, and 30 for depths z1, z2, and zN, respectively. Migration operator M(xcfx89,{right arrow over (xcex7)},z) is then 2-D Fourier transformed for every frequency xcfx89 and every depth z into the frequency-wavenumber (xcfx89xe2x88x92{right arrow over (k)}) domain resulting in transformed migration operator {tilde over (M)}(xcfx89,{right arrow over (k)},z). In FIG. 2, the transformed migration operator is represented by cubes 32, 34, and 36 for depths z1, z2, and zN, respectively. The transformed migration operator {tilde over (M)}(xcfx89,{right arrow over (k)},z) for each depth is then multiplied by the transformed data {tilde over (D)}(xcfx89,{right arrow over (k)}). For each depth z1, z2, . . . zN, the migrated data are then summed over frequency to form an image {tilde over (D)}M({right arrow over (k)},z) in the wavenumber domain. In FIG. 2, this is represented by images 38, 40, and 42 for depths z1, z2, and zN, respectively. Finally, the wavenumber domain migrated images {tilde over (D)}M({right arrow over (k)},z) are 2-D Fourier transformed from the wavenumber domain back to the space domain resulting in migrated data DM({right arrow over (x)},z), which is represented in FIG. 2 by migrated data cube 44.
The method described by Etgen is substantially faster than Kirchhoff prestack migration according to equation (1), but is dominated by the costs of computing and transforming the migration operator repeatedly for each depth. As further described below, in some embodiments of the present invention the unmigrated data traces are transformed to the xcfx89xe2x88x92{right arrow over (k)} domain in the same manner as proposed by Etgen; however, the present invention computes and applies the migration operator in a unique way that requires significantly less computer resources than the Etgen method.
From the foregoing, it can be seen that there is a need for a method for economically achieving accurate prestack migration of seismic data. Such a method should be capable of accurately preserving the seismic amplitudes of the unmigrated data. Such a method should also be capable of imaging reflectors having very steep dips, including reflectors having dips above 90 degrees dip angle such as may be found, for example, adjacent to salt structures in the Gulf of Mexico. The present invention satisfies this need.
In one embodiment, the present invention is a method for prestack time migration of common-offset seismic data traces which comprises the steps of (1) selecting a seismic velocity function for use in the prestack time migration, (2) transforming the common-offset seismic data traces from the space-time ({right arrow over ("xgr")}xe2x88x92t) domain to a preselected alternate data domain, (3) using the seismic velocity function to compute a time function in the {right arrow over (p)}xe2x88x92z domain for accurately migrating the seismic data traces, (4) computing a weighting function in the {right arrow over (p)}xe2x88x92z domain, which weighting function substantially preserves the seismic amplitudes of the seismic data traces, (5) using the time function to phase shift the seismic data traces in the alternate data domain, (6) multiplying the amplitudes of the phase-shifted seismic data traces by the weighting function, (7) summing the phase-shifted and amplitude-weighted data over an appropriately chosen variable so as to form an image in said alternate data domain, and (8) performing a 2-D transformation of the image from the alternate data domain back to the space ({right arrow over (x)}xe2x88x92z) domain.
The preselected alternate data domain may be the xcfx89xe2x88x92{right arrow over (k)} domain, the {right arrow over (k)}xe2x88x92t domain, or the xcfx84xe2x88x92{right arrow over (p)} domain.
The inventive method also permits imaging of seismic turning ray energy for non-zero source-receiver offsets.